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The development and application of a multidimensional numerical simulation code for investigating near- 
field plasma processes in magnetic nozzles are presented. The code calculates the time-dependent evolution 
of all three spatial components of both the magnetic field and velocity in a plasma flow, and includes physical 
models of relevant transport phenomena. It has been applied to an investigation of the behavior of plasma flows 
found in high-power thrusters, employing a realistic magnetic nozzle configuration. Simulation of a channel- 
flow case where the flow was super-Alfvenic has demonstrated that such a flow produces adequate ‘back-emf 9 
to significantly alter the shape of the total magnetic field, preventing the flow from curving back to the magnetic 
field coil in the near-field region. Results from this simulation can be insightful in predicting far-field behavior 
and can be used as a set of self-consistent boundary conditions for far-field simulations. Future investigations 
will focus on cases where the inlet flow is sub-Alfvenic and where the flow is allowed to freely expand in the 
radial direction once it is downstream of the coil. 


Nomenclature 


B 

magnetic induction, T 

P 

pressure, Pa 

Sm 

Maxwell stress tensor 

q 

dissipative power per area, W/m' 

E 

electric field strength, V/m 

T 

temperature, K 

E f 

electric field in the plasma frame, V/m 

u 

fluid velocity, m/s 

£ 

energy density of the plasma, J/m 3 

rj 

resistivity, H.m 

k 

thermal conductivity, W/m.K 

po 

permeability of free space, H/m 

j 

current density, A/m 2 

p 

mass density, kg/m 3 


I. Introduction 

M agnetic nozzles are used in many plasma thrusters to convert random thermal energy into directed kinetic 
energy and to control the direction of the exhaust. However, due to the solenoidal nature of the magnetic field 
it is essential to ensure that the plasma detaches from the field, making this a problem of active research. 1 ’ 2 While 
other investigations have looked into far-field phenomena, like detachment, this paper discusses numerical simulations 
of the near-field behavior of plasma jets in magnetic nozzles, where the applied field gradients are greatest. The 
development of a multidimensional numerical tool that can be used to investigate plasma detachment mechanisms 
was introduced in Ref. [3]. This code solves the conservation form of the magnetohydrodynamic (MHD) equations, 
accounting for thermal nonequilibrium between electrons and the rest of the plasma and including appropriate models 
for ionization, equation of state and various transport phenomena. This is an extension of previous simulations of 
self-field magnetoplasmadynamic thrusters (MPDT) 4 ’ 5 and a lithium Lorentz force accelerator. 6 

The outline for the rest of this paper is as follows. The governing equations for the nozzle simulations are described 
in Sect. II. The geometry of the model, including the description of the magnetic nozzle, and the boundary conditions 
used are discussed in Sect. III. We present the results and a discussion of a plasma simulation test case involving a 
super-Alfvenic channel flow in solenoidal magnetic fields in Sect. IV, and discuss the next steps of this effort in Sect. 
V. 
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II. MHD Equations for Simulations 


The governing equations for an MHD flow problem are represented in conservation relations for mass, momentum, 
magnetic flux, and energy. These can be compactly written in the form, 
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puu + p- B m 
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£ 


_ (£ + p) u - B m • u _ 


The right hand side of the equation, D, contains dissipative effects, which are also written as the divergence of fluxes. 6 
The compact vector-tensor equations are expanded into their individual vector components and rewritten in a finite- 
volume form so that they can be solved. The expanded equation set, assuming no azimuthal variation (< d/80 = 0), 
is: 
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Qr = (K b 6 ~ E e B d /Mo + ( k rr dT e /dr ) + (fc rz dT e /dz) , 
q z = ( K B e ~ K B e ) /Mo + ( k zr dT e /dr ) + (fc zz dT e /dz) . 


The thermal nonequilibrium between the electrons and the heavy species (ions of varying levels of ionization and 
neutrals) is accounted by maintaining a separate energy equation for the electrons. In the presence of an applied 
magnetic field, the components B r and B z are the sum of the steady-state applied field and an induced field that is 
initially time-varying but eventually converges to a steady-state value. The set of MHD equations is closed by the 
generalized Ohm’s law, 


E' = E + (u x B) = 77J + 


(j X B - Vp e ) 

en e 


( 3 ) 


and Ampere’s law. The models for transport coefficients, equation-of-state, ionization and thermal nonequilibrium are 
described in Ref. [6]. 
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III. Simulated Geometry and Boundary Conditions 


m. A. Physical Domain and Magnetostatic Model 

A separate magnetostatic solver 7 is used to create an applied magnetic field topology. The magnetic coils in real 
experiments have tens of thousands of loops, while the coils in magnetostatic simulation had only a handful of loops. 
Consequently, the primary goal of our magnetostatic modeling effort was to recreate a reasonable and realistic topol- 
ogy. In Fig. la we show the nozzle topology for the experiment described by Kuriki and Okada 8 generated using the 
magnetostatic solver. 



Figure 1. a) Schematic of the experimental setup (based upon Ref. [8]) and the domain that was modeled in this paper. Magnetostatic simulation results 
showing the magnetic flux lines in the r — z plane of the nozzle are presented, with these data used as the applied field background in this simulation; b) 
Magnitude of the applied magnetic field produced by the coil. 


The computationally-generated field was compared to published values from the experimental setup. Despite our 
simplifying the problem by using only a few magnet coil loops, the magnetostatic solver produced a magnetic field 
topology possessing field values that were within a few percent of those measured. 8 The magnetostatic solver 7 yields 
the relevant variables (B r , B z ) on an orthogonal (r, z) grid. An interpolation script was used to transfer the magnetic 
field data onto the non-orthogonal grid used by the MHD code. Though it is not our goal to simulate the experiment 
in Ref. [8], our simulation uses the magnetic field topology used in that experiment because of its simplicity and 
appropriateness for our investigation. 

The domain of the simulation is shown in Fig. la as the dotted rectangle. Here, the bottom surface of the coil is at 
a distance R = 5.0 cm from the centerline and has a 5.0 cm by 5.0 cm rectangular cross-section, has a thickness of 5.0 
cm as well. The computational domain extends 5.0 cm upstream and 10.0 cm downstream from the center of the coil. 
The domain is bounded at the top (i2 ma x = 4.98 cm) by an insulator, resulting in a restricted channel-flow simulation. 

m.B. Boundary Conditions 

Near-field investigations in applied solenoidal fields pose special challenges for numerical simulations owing to the 
strong curvature of the magnetic field, large gradients in the flow parameters, and interactions with material boundaries. 
We proceed with a discussion of the boundary conditions used in this paper. 

III.B. 1 . Flow Properties 

• Inlet - At the inlet (z = 0), a specified mass flow rate of the propellant enters at a specified temperature at sonic 
conditions. The velocity vector is set to be parallel to the magnetic field vector at the inlet, treating the incoming 
plasma as fully magnetized. 

• Exit - At the axial outflow boundary, the axial gradient of the flow quantities is set to be same as the axial 
gradient just inside the exit plane. 

• Solid Walls - In the radial direction, at # max , the flow is obstructed by a solid surface. Here, we enforce 
n ■ u = 0. This surface is treated as a thermal insulator. Consequently thermal conduction is set to zero. 
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• Centerline - At the axis of symmetry (r = 0), there are no radial convective or conductive fluxes. Moreover, 
there are no radial gradients here. 


III.B.2. Field Properties 

• Inlet and Exit - In this simulation, the computational domain is not large enough to assume that the magnetic 
field is zero at the boundaries. In fact, such an assumption is not justified given that one of the primary objectives 
of this investigation is to study the interplay between the plasma and the magnetic nozzle. Therefore, the 
magnetic field in the ghost cells adjacent to the inlet and the exit consists of both the applied field of the 
magnetic nozzle and the induced field produced by internal plasma currents. Contributions from the latter 
are self-consistently calculated for each ghost cell (at the general location (r g , z g )) using the Biot-Savart law, 
summing the contributions from all the currents inside the domain, 9 

(4) 

While more efficient ways for implementing this can be found in Ref. [10], we employ an explicit numerical 
evaluation of Eq. (4). 

• Solid Walls - At the insulated boundary at i? max , the magnetic field is assumed to diffuse instantaneously into 
the wall (i.e., rj — ► oo) implying that the surface current is zero. 

• Centerline - At the axis of symmetry j r = jo = 0, and the inductive component of the electric field is zero 
because u r = uq = B r = Bq = 0. The axial current density, j z can be obtained from the magnetic field at the 
point adjacent to r = 0 using a Taylor series approximation, which yields 


Jz\ r= Q ~ 


4 Bq\ Ar 
1 2 

/i 0 Ar 


( 5 ) 


IV. Results and Analysis 


Theories of plasma detachment from magnetic fields 11, 12 identify the ratio of the flow speed to the AHVen speed, 
va = B/^/fji 0 p 9 as a key parameter in this process. This ratio, also known as the AHVen Mach number M A , cor- 
responds to the square root of the ratio of the kinetic energy of the plasma to the magnetic field energy which, in 
equi librium , scales with the square root of the plasma beta, 


M a 
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V A 
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B/ yj pop 


f™ 2 / 2 ~ / P - Jq 

B 2 j (2/i 0 ) — Y B 2 / (2/i 0 ) VP ' 


( 6 ) 


We chose to examine a channel flow case where the 
flow speed exceeds the AHVen speed to analyze the how 
the flow and the field shape each other. Therefore we set 
uf v A — M a = 2.0, based on a magnetic field value of 
301 G. This corresponds to the value of the initial ap- 
plied field under the center of the coil (z = 5 cm) taken 
at the radial location r = 3.4 cm, which corresponds to 
the radius where half the mass in a uniform flow would 
be below this radius and half above. The mass flow rate 
was fixed at m = 6.0 g/s, and the temperature at the inlet 
was set at T e = Th = 1.5 eV. These parameters deter- 
mine the density and the speed of the incoming plasma 
flow. The plasma conditions chosen are representative of 
those found in propulsive flow exhausts. For a domain 
that extends 1 5 cm in the axial dimension, a plasma flow 
of u z ~ 3000 m/s takes 0( 100)ps to reach establish 
equilibrium. 

Fig. 2 shows the ratio of the flow speed to the AllVen 
speed at the end of the calculation. Since the flow is set to 



M a = 1.7 6.5 11.0 

Figure 2. Profile of the Alfven Mach number, Ma, obtained for a converged, 
steady-state solution. 

be super-Alfvenic at the region of strongest magnetic field, 
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directly below the coil, Ma is even higher in the regions where the magnetic field is weaker. Under these conditions 
the magnetic Reynolds number Re m ~ 0(1), indicating that convective effects are significant, and that the ’back 
emf will play an important role in this problem. This is consistent with the results of this simulation. With Ma > 1, 
the induced field shields out the applied field. This is primarily accomplished by large azimuthal currents that are 
induced in the plasma. As shown in Fig. 3, jo is large below the comers of the coil where the gradients in the applied 
magnetic field are the strongest. The large azimuthal current in a thin layer near the inlet is an artifact of the boundary 
counditions, due to which the incoming flow encounters the magnetic field and begins to swirl. 

Because Rem ~ 0(1), dissipative effects are also 
significant. Consequently, the magnetic field is not com- 
pletely shielded from the interior of the plasma, with 
some of the field diffusing inwards. In Fig. 4, flux tubes 
of the initial (applied) and the final (applied + induced) 
magnetic field are presented, showing how the plasma 
acts to shape the total magnetic field. Due to the large 
azimuthal currents induced in opposite directions at the 
edges of the coil, shown in Fig. 3, the magnitude of the 
axial magnetic field between the bounds of the coil is 
very low, and the field downstream of the coil is signifi- 
cantly altered by the plasma. 

The final field is stretched axially at the exit, as 
shown in Fig. 4, allowing the plasma to exit the domain 3e =- 21 - 10 0 10 2 ikA/m 

without closing back to the coil in the near-field region. 

The angle made by the initial (applied) and final (applied Figure 3. Azimuthal plasma currents, j e , calculated within the domain. 

+ induced) fields with the z axis can be a useful metric to understand the effect of the plasma flow in shaping the mag- 
netic field. As shown in Fig. 5, the axial flux of the plasma out of the domain is not seriously impeded by the applied 
magnetic field. In fact, in this case the field and plasma align themselves with each other quite well, with both vectors 
pointing primarily in the axial direction at the exit. 




Figure 4. Magnetic flux tubes of a) the initial (applied) magnetic field and b) the final (total) magnetic field, showing the effect of a super-Alfvenic plasma on 
the field distribution. 


V. Concluding Remarks 

We have presented the development and application of a code that calculates the time-dependent evolution of 
all three spatial components of magnetic field and velocity of plasma flows in magnetic nozzles. We have used it 
to investigate the behavior of a plasma flow, with plasma parameters typically found in high-power thrusters, in a 
magnetic nozzle configuration used in an experiment. We simulated a channel-flow case where the flow was super- 
Alfvenic, and demonstrated that such a flow produces adequate ‘back-emff to shape the total field and prevent the 
flow from curving back to the magnetic coil in the near-field region. We showed that this is accomplished through the 
production of large azimuthal currents, especially below the coil comers where the applied field gradients are largest. 
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Figure 5. The angle made by the initial (applied) and final (applied + induced) fields with the z axis at the exit plane. 


In the future, we intend to simulate cases where the inlet flow is sub-AHVenic to compare with the results of the present 
simulation, and also allow unimpeded radial expansion of the plasma jet downstream of the coil. 

Given the limitations of the continuum model used in this work, we have restricted our investigation to the near- 
field region of the magnetic nozzle. Other investigations, such as that found in Ref. [1], discuss the far-field region. 
Nevertheless, observations of plasma behavior in the near-field of the magnetic nozzle have proved to be informative 
and the results from this simulation can be extrapolated to estimate far-field behavior and can yield self-consistent 
boundary conditions for far-field simulations. 
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